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Abstract 


-- lt  > 


This  paper  describes  a method  of  generating  gamma  variates  that 

otk**  r 


appears  to  be  less  costly  than'a  recently  suggested  method  In  [3].  For 
large  shape  parameter  a' the  cost  of  computation  is  proportional  to 

- ^ * pf-4 

whereas  the  „ method r\ n [3]  Is  proportional  to  a.  Experimentation^ in  [2] 

<o  indicates  that  for  small  a the  method  suggested  here  also  dominates 

methods  recently  suggested /fn"  tU.^albeit  those  methods  dominate  for  large 

* V ; * 

a.  The  method  suggested  here  uses  the  rejection  technique. 


1 . Introduction* 

This  paper  describes  a new  technique  (method  1)  for  sampling  from 
the  gamma  distribution  on  a digital  computer  and  compares  It  with  an 
alternative  technique  (method  2)  that  Wallace  has  suggested  In  [3].  A 
gamma  variate  X has  the  probability  density  function**  (p.d.f.) 
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(1)  fx(x) 


xa_1  e'x/r(a) 
0 


0 < X < a,  a > 0 
elsewhere. 


Both  methods  use  the  rejection  method  and  apply  for  a > 1. 

2.  Rejection  Method 

Let  X be  a nonnegative  valued  continuous  random  variable  with 
bounded  p.d.f.  representable  in  the  form 


C(a,  e)a(a,  e)g(x,  a,  g)h(x,  a,  B) 

0 5.  x 1.  " 

(2) 

f v ( X j a)  = J 

0 

elsewhere 

0 <_  h(x,  a,  B ) > 

00 

f h(x,  a,  B)dx  = 1 , 
0 

0 < g(x,  a,  b) 

< a(a,  B)  > l/g(x,  a,  B), 

l/c(a,  B)  = a(a 

CO 

, B)  / g(x,  a,  B)h(x,  a,  B)dx. 

0 


Let  X'  denote  a random  variable  with  p.d.f.  h and  let  U be  a uniform 
deviate  on  (0,  1).  If  U£a(o,  e)g(X',  a,  6)  then  X'  has  the  p.d.f. 


+ I am  grateful  to  Mr.  Hunter  McDaniel  for  programming  methods  1 
and  2 in  PL/1. 

Here  we  assume  a unit  scale  parameter  without  loss  of  generality. 
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fx  In  (2).  This  result  follows  from 


f X , (x I U _<  a(a,  B)g(x,  a,  B))  = 
(3) 


pr[U  < a(o,  B)g(x,  a,  e) | X ' = x]h(x,  a,  b) 
pr[U  < a(a,  B)g(x,  ot,  B)] 


“ fx(x,  tt ) • 


Since 


(4)  pr[u  < a(a,  B)g(x,  a,  B)]  = l/c(a,  B) 


c(a,  b)  denotes  the  mean  number  of  trials  to  obtain  an  X from  (2).  For 
a given  X'  from  a specified  h we  want  the  probability  of  success  to  be 
as  close  to  unity  as  possible.  This  feature  requires 


(5)  l/a*(a,  B)  = max  g(x,  a,  s). 


For  any  X'  we  want  (4)  to  be  as  large  as  possible,  which  implies. 


c(a,  s*)  = min  c(a,  b)  = min[l/a*(a,  e)Eg(X 1 , a,  b)] 

B 6 

= min[max  g(x,  a,  B)/Eg(X',  a,  b)] 
B x 

00 

Eg( X ' , a,  B)  = / g(x,  a,  B)h(x,  a,  B)dx. 

0 


The  distinction  between  methods  1 and  2 lies  in  the  choice  of  h. 

Table  1 shows  relevant  quantities  for  each  proposal.  To  make  an  appropriate 


comparison  between  methods  we  need  to  consider  the  mean  number  of  trials 

c . (a,  8*)  for  each  and  the  mean  number  of  required  random  numbers. 

J 


Table  1 

Gamma  Generation*  for  a > 1 


Method 

i 

hj(x,  a,  8) 

g^x,  a,  e) 

a.j*(a,  8) 

B1* 

C.j (a,  B*) 

1 

B-V*/e 

xa-le(l/B-l)x 

(e/a)”-1 

a 

a“e1-a/r(a) 

2 

xY-le-x[(l_6)Y+8x] 

XY< 

yO-b) 

b(1-y')’ 

y' 

Y* 

r(Y)Y1-Y,/r(Y) 

r(Y+l) 

Y=<a> 

(1-B)y+BX 

y'sa-y 

l-Y* 

y(1-b)y' 

* <e>  denotes  the  largest  integer  in  e. 

3.  Method  1 

Conceptually,  method  1 Implies  4 steps: 

1.  Generate  an  exponential  deviate + X'. 

2.  Generate  a uniform  deviate  U. 

3.  If  U ^ (X'/eX  +1)a-1  then  X = aX'  has  the  p.d.f.  in  (2). 

4.  Otherwise,  return  to  step  1. 

If  we  use  the  inverse  transform  method  to  generate  X'  then  each  trial  requires 
2 random  numbers.  Therefore,  the  mean  number  of  random  numbers  needed  to 
generate  X from  (2)  is  2aa/r(a)ea"1 . For  large  a this  quantity  is  approximately 
e(2a/w)1^2,  an  appealing  result.  Notice  that  for  large  Integral  a,  using  method  1 


t 


An  exponential  deviate  has  unit  mean. 


requires  fewer  random  numbers  than  the  conventional  method  which  uses 

b 

(7)  X = -Inn  (Uj, 

i-1  1 

U,,  ....  U being  a sequence  of  independent  uniform  deviates.  For  small  integral 
a one  can  show  that  (7)  is  superior.  Our  experiments  indicate  that  method  one 
prevails  for  nonintegral  a < 7 and  all  a > 7. 

4.  Method  2 

For  integral  a method  2 uses  (7).  For  nonintegral  a the  6 steps  are: 

1.  Generate  a uniform  deviate  U. 

2.  If  l)  <_  1 -a+<a>  generate  X1  from  (7)  using  b = <a>. 

3.  Otherwise,  generate  X'  from  (7)  using  b = <a+l>. 

4.  Generate  a uniform  deviate  U. 

5.  If  U < (X'/y)y  /(l-y'+X'/y)  then  X'  has  the  p.d.f.  (2). 

6.  Otherwise,  go  to  step  1. 

These  steps  require  a+2  random  numbers  on  average  per  trial.  Therefore, 
for  nonintegral  a,  method  2 uses  (a+2)r(Y)Y^  y /r(a)  random  numbers  on  average. 
This  quantity  is  approximately  a+2  for  a>5. 

5.  Comparison  of  Methods 

PL/1  programs  were  prepared  using  algorithm  G1  for  method  1 and 
using  the  steps  given  in  [3]  for  method  2. 

Algorithm  G1 


Given:  a 

1 . a'  *■  a-1 . 

2.  Generate  a uniform  deviate  U. 
(continued) 
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3.  V «-  -In  U. 

4.  Generate  a uniform  deviate  U. 

5.  W <-  -In  U. 

6.  If  W >_  a ' ( V - In  V-l ),  X aV  and  return  with  X. 

7.  Otherwise,  go  to  2. 

Table  2 displays  the  results  for  generation  of  10,000  gamma  variates  for  each 
selected  value  of+  a. 

Table  2 

Comparison  of  Methods 


mean 

(in 

CPU  time 
psec. ) 

a 

Method  1 

Method  2 

Ratio 

1.25 

723 

1093 

.661 

2.25 

988 

1300 

.760 

3.25 

1193 

1542 

.774 

4.25 

1352 

1787 

.757 

5.25 

1541 

2039 

.756 

Based  on  these  results  we  computed  expressions  for  T^,  the  mean 
CPU  time  for  method  i,  as  a function  of  a.  These  expressions  are  in  microseconds. 

T1  - 140  + 624aa/r(a)ea'1 

(8) 

T2  = -88+(752+254a)r(y)y1'a+Y/r(a). 


The  programs  were  run  on  the  IBM  360/75  computer  at  the  University  of  North 
Carolina  Computer  Center  at  Chapel  Hill  as  single  stream  inputs.  This 
procedure  minimized  the  error  due  to  monitoring  in  a multiprogram  mode. 
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For  large  a T-|/T2  *v  624/254/x  = 2.46//x.  For  example,  a = 30  gives 
T i /T2  ^ 0-45  and  a = 50,  T-j/T2  ^ 0.35. 

One  modification  to  method  1 makes  it  at  least  as  good  as  method  2 for 
all  integral  a,  while  preserving  its  superiority  for  nonintegral  a.  Exper- 
imentation with  method  1 revealed  that  it  is  superior  to  method  2 for  all 
a > 7.  Addition  of  the  statement: 


a 

0.  If  a £ 7 and  <a>  = a,  return  with  X = -ln(n  U.) 

i=l  1 

prior  to  statement  1 in  algorithm  G1  modifies  the  flow  appropriately. 

5.  New  Prospects 

Upon  conclusion  of  the  work  presented  here  the  writer  learned  of  research 
by  Dieter  and  Ahrens  in  [1]  on  gamma  generation  1)  using  a truncated  noncentral 
Cauchy  distribution  for  h and  2)  exploiting  the  relationship  between 
the  gamma  and  normal  distributions  for  large  a.  The  most  notable  feature 
of  their  work  is  that  computation  time  goes  to  a fixed  limit  as  a increases. 
Although  this  property  makes  the  Dieter  and  Ahrens  procedures  more  attractive 
for  large  a,  Robinson  and  Lewis  [2]  have  recently  prepared  a gamma  generation 
program  in  which  a variant  of  algorithm  G1  dominates  all  competitors  for 
1.2  <_a  £2.9.  Since  this  is  a commonly  encountered  range  in  practice  the 
significance  of  method  1 remains. 

Since  the  work  in  [2]  generates  exponential  variates  by  a more  efficient 
method  than  inverse  transformation  does,  it  is  not  presently  clear  to  the 
writer  what  the  range  of  superiority  would  be  using  algorithm  Gl.  This  issue 
is  a legitimate  concern  since  simulation  languages  such  as  SIMSCRIPT  and 
SIMPL/1  use  the  inverse  approach. 
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